package scu.maqiang.homogenization;

import scu.maqiang.fes.BVPType;
import scu.maqiang.fes.FES2Q41;
import scu.maqiang.mesh.Mesh2Q4;
import scu.maqiang.mesh.ScalarFunc;
import scu.maqiang.numeric.Direct;
import scu.maqiang.numeric.NewIterSSolver;
import scu.maqiang.numeric.SRMatrix;

public class Possion2DSOTSHole {

	public static void main(String[] args) {
		Mesh2Q4 mesh = new Mesh2Q4().square2D(256, 256);
		int n = 16;
		double eps = 1.0/n;
		ScalarFunc chFun = (xy, label, param) -> {
			double xCell = xy[0] * n - (int)(xy[0] * n);
			double yCell = xy[1] * n - (int)(xy[1] * n);
			return (Math.abs(xCell- 0.5) < 0.25 && Math.abs(yCell - 0.5) < 0.25)? 1.0:-1.0;	
		};
		
		mesh.removeElements(chFun, null, 8);
		FES2Q41 fs = new FES2Q41(mesh);
		SRMatrix A = new SRMatrix(fs.GetNdof());
		double[] RHS = new double[fs.GetNdof()];
		fs.assembleStiff(new double[] {1.0}, BVPType.COMMON, A);
		fs.assembleSource(new double[] {1.0}, BVPType.COMMON, RHS);
		fs.applyBC_MBN(A, Direct.All, 1, 2, 3, 4);
		fs.applyBC_MBN(RHS, Direct.All, 0.0, 1, 2, 3, 4);
		NewIterSSolver solver = new NewIterSSolver(A);
		double[] x0 = new double[fs.GetNdof()];
		solver.PCGSSOR(RHS, x0, 1.5, 1);
		mesh.toTecplot("PossionHole.dat", x0);
	}

}
